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Stochastic differential equations (sdes) play an important role in physics but existing numerical 
methods for solving such equations are of low accuracy and poor stability. A general strategy for 
developing accurate and efficient schemes for solving stochastic equations in outlined here. High 
order numerical methods are developed for integration of stochastic differential equations with strong 
solutions. We demonstrate the accuracy of the resulting integration schemes by computing the errors 
in approximate solutions for sdes which have known exact solutions. 
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Stochastic differential equations (sdes) have a long his- 
tory in physics[l| and play an important role in many 
other areas of science, engineering and finance^, 0, 0J- 
Recently a number of computational techniques have 
been developed in which high dimensional determinis- 
tic equations are decomposed into lower dimensional 
stochastic equations. Gisin and Percival4], for exam- 
ple, reduced a deterministic master equation for the den- 
sity matrix into stochastic equations for a wavefunction. 
Similar approaches are being used to solve the quan- 
tum many-body problem for bosons 0, fermions0 and 
vibrations 0. These latter methods give rise to large sets 
of coupled sdes which require fast and efficient numerical 
integration schemes. Unfortunately, and in spite of their 
widespread use, the available numerical techniques^ for 
solving such equations are far less accurate than compa- 
rable methods for solution of ordinary differential equa- 
tions (odes). 

In this manuscript we show how classical methods for 
solving odes, such as Runge-Kutta, can be adapted for 
the solution of a class of sdes which should include many 
of the equations which arise in physical problems. 

Consider a finite set of sdes, 
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represented in Ito jUSOl form, where j — 1, . . . , n. Here 
X t = (Xt,...,Xt) and the dW k are independent and 
normally distributed stochastic differentials with zero 
mean and variance dt (i.e. sampled N(0,dt)). The 
stochastic variables W 4 are Wiener processes. Now as- 
sume that the coefficients a 3 and b 3 k have regularity prop- 
erties which guarantee strong solutions, i.e. that X\ are 
some fixed functions of the Wiener processes, and that 
they are differentiable to high order. [Sufficient con- 
ditions for strong solutions are discussed in Ref. [jj.] 
We may then view the solutions of Q as functions 
X{ = Xj(t, W t \..., W™) of time and the Wiener pro- 
cesses. The solutions can therefore be expanded in Taylor 
series. Keeping terms of order dt or less then gives 
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In a mean square sense the product of differentials 
dW^dWj is equivalent to Skddt in the It60,0,lj| formu- 
lation of stochastic calculus. Making this replacement 
then yields 
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which when compared to (JIJ allows us to identify the first 
derivatives 
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Now that these first order derivatives are expressed in 
terms of a? and bl, higher order derivatives can be com- 
puted. Thus a Taylor expansion of the solutions 
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can be obtained for finite displacements At and AW k . 
This Taylor expansion can then be employed to develop 
Runge-Kutta algorithms and other integration schemes. 

We illustrate the use of this approach by developing 
a Runge-Kutta method for sdes which is closely related 
to the classical Runge-Kutta scheme for odes. For given 
displacements At and AW k define 
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and consider the following four stage approximation 
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where t, is the initial time and tj+i = ij + At. Taylor 
expansion of this scheme shows that X ti+1 differs from 
the exact solution by terms of order higher than At 2 (i.e. 
terms of higher order than At 2 , At(AW t fe ) 2 , (AW?) 4 , 
(AW t k ) 2 (AW[) 2 , and (AW t k ) 2 AW l t AW/) . Thus, this 
stochastic Runge-Kutta algorithm plays a role very sim- 
ilar to its classical counterpart except that its order is 
reduced from four to two. Generalizations to higher or- 
der Runge-Kutta schemes are straightforward, and we 
will employ one such scheme in example calculations, but 
details will not be presented here. 

While this approach is not completely general, since it 
will fail for sdes with weak solutions or non-differentiable 
a? and b J k , it should be applicable to a wide range of prob- 
lems. It can for example be used to solve every one of the 
equations with known solutions tabulated in section 4.4 
of Ref . 3] . To illustrate the accuracy of the method and 
its improvement over other known techniques for solv- 
ing sdes we now consider a number of these examples. 
We compare known exact solutions with numerical so- 
lutions obtained using the Euler-Maruyama scheme^, 
a derivative free version of the Milstein scheme due to 
Kloeden and PlatenQ, the classical Runge-Kutta scheme 
(JSJ, and another Runge-Kutta scheme obtained in the 
manner outlined above from an eighth order twelve step 
method for odes due to Hairer and Wanner [Io| (this re- 
produces the stochastic Taylor expansion up to and in- 
cluding terms of order At 4 ). Stochastic differentials were 
sampled using the routines gasdev and ran2|ll|. 

As a first test of these methods consider an au- 
tonomous nonlinear scalar equation 

dX t = (1 + X t )(l + X 2 )dt + (1 + X 2 )dW t (9) 

with just one Wiener process. In this example and in all 
subsequent examples we assume all Wiener processes are 
initially zero. The exact solution to this equation isQ 

X t = tan(t + W t + arctan(Jfo)) (10) 

as can be readily verified using ItoQ, 0, calculus. In 
Fig. 1 we plot the error log 10 \X t — x^ pproxima e | vs 



Figure 1: log 10 \Xt ~ X? pproxlmate \ vs time t for Eq. © 
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Figure 2: log 10 \X t - x? pproxlmate \ vs time t for Eq. JTTJ 
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time computed with a time step of 2.5 x 10 for a sin- 
gle stochastic trajectory with initial condition X a = 1 
for the four different approximation schemes. The Mil- 
stein scheme (long-dashed curve) shows some improve- 
ment over the primitive Euler-Maruyama method (solid 
curve) but the order two Runge-Kutta scheme (short- 
dashed curve) and order four Runge-Kutta scheme (dot- 
ted curve) perform very much better. 

The second example equation, also from Ref. 0, is 
an autonomous linear scalar equation in two Wiener pro- 
cesses 

dX t = a X t dt + b x X t dW} + b 2 X t dW 2 (11) 
which has an exact solution 

X t = X exp{[a - i(6 2 + 6 2 )]t + b x W} + b 2 W 2 }. (12) 

The logarithm base ten of the error for the different 
schemes, calculated for initial condition Xq = 1 and time 
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step .01, is plotted in Fig. 2. Here the Milstein scheme 
(long-dashed curve) performs no better than the Euler- 
Maruyama method (solid curve) but again the order two 
Runge-Kutta scheme (short-dashed curve) and order four 
Runge-Kutta scheme (dotted curve) show greatly im- 
proved accuracy. [Note that the apparent improvement 
in performance of all schemes at long time is a result of 
the fact that the solution decays to zero.] 

Example 3 is a set of two coupled linear autonomous 
sdes 
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with three Wiener processes. Here the solutions are 
X\ = exp{-2t + Wl ~ W? }cosFT t 3 



X 2 = exp{-2t 
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Numerical solutions were calculated with a time step of 
.01 and errors in X\ are represented in Fig. 3. The order 
two Runge-Kutta scheme (long-dashed curve) and order 
four Runge-Kutta scheme (short-dashed curve) show im- 
provement over the Milstein scheme (solid curve). Simi- 
lar results were obtained for X 2 . 

The examples we have considered so far have not had 
explicitly time dependent a J and b J k . Example 4 is a 
scalar non-autonomous sde 
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with known solution^ 



X t = 



' ' ' X + ^(l + t) 2 (W t + t-t ). (16) 
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Numerical solutions were calculated using the order two 
Runge-Kutta scheme and a time step of .001, to = 
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and Xq — 1. The error is represented in Fig. 4. As in 
previous examples a high accuracy is achieved in spite 
of the rapid growth of the solution. The comparative 
smoothness of the error curve reflects the fact the the 
deterministic part of the solution dominates. 

We now consider an example for which an exact solu- 
tion is known but which is expressed in terms a stochastic 
integral. Consider the stochastic Ginzburg-Landau equa- 
tion 



dX t = [~Xf + (a+ ^a 2 )X t ]dt + crX t dWt 



with solution^ 
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We chose a = .01, a = 4, X = 1 and dt = 5 x 10~ 6 . The 
stochastic integral was computed using a Riemann sum 
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Figure 6: log 10 \X t - X ? pproxlmate \ vs time t for Eq. QSJ 
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with exact solution 3] 

X t = arcsinh [er at sinh X + e~ at J e as dw}j . (20) 

We set a = .02, b = 1, X = 1 and dt = 1 X 10~ 5 . The 
stochastic integral in the exact solution was calculated 
using the It6[l|,|2j,|3j integral formula with the same time 
step. The error in the solution calculated with the order 
two Runge-Kutta scheme is plotted in Fig. 6. As in all 
previous cases considered the accuracy is very good. 

Thus, the approach to solving sdes advocated here 
works very well for the wide range of examples we 
have considered. The order 4 Runge-Kutta method is 
clearly much more accurate than the order 2 Runge- 
Kutta scheme. It also has an embedded lower order 
Runge-Kutta scheme which can be employed to obtain 
an error estimate suitable for stepsize control [loj. Hence 
is should be possible to use variable stepsizes to ensure 
the accuracy of the solution. This sort of implementa- 
tion is essential for solving equations which do not have 
known exact solutions. The only subtlety in developing 
such a method is ensuring that the correct Wiener path 
is maintained even when a step must be rejected. This is 
achieved by dividing the rejected differentials dt and 
dWf in two segments; dt/2 and dWf /2 — y followed by 
dt/2 and dWf /2 + y where y is sampled N(0, dt/2). To 
illustrate the accuracy of the resulting variable stepsize 
algorithm we solve the Gisin-PercivalQ stochastic wave 
equation for the nonlinear absorber (Eq. 4.2 of Ref. 0) 

d\ip) = .l(a f - a)\tp)dt + (2a^a 2 - a t2 a 2 -a^ aS)\i/f)dt 
+ V2(a 2 -H?)\il))dW t (21) 



with the same time step. Error in the solution calculated 
with the order two Runge-Kutta scheme is plotted in Fig. 
5. Good accuracy is again obtained. 

Finally, we consider an example in which the exact 
solution is expressed in terms of a ItopJ 0, stochastic 
integral. Consider the sde 

dX t = -t&nhX t (a+^b 2 sech 2 X t )dt + bsechX t dW t (19) 



with initial state |^(0) >= |0). In Fig. 7 we plot the 
error in mean occupation number n t — M[{i/>\a* a\i/>)] vs 
time (Fig. 5 of Ref. 0) where M[-\ denotes an average 
over stochastic realisations. 1000, 10000, and 20000 tra- 
jectories were used to calculate the solid curve, dashed 
curve and dotted curve, respectively. Convergence to the 
exact result is good. 
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